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Abstract.  The  generation  of  vortices  for  the  flow  past  a  flat  plate  is  calculated  by  the  DSMC  method.  The  new 
intermolecular  collision  scheme  and  the  conventional  one  are  compared  with  respect  to  their  allowable  cell  dimension.  In 
the  condition  of  Re= 50,  when  using  the  conventional  scheme,  the  vortex  behind  an  inclined  flat  plate  becomes  weak  with 
increasing  the  cell  dimension  beyond  the  mean  free  path,  and  finally  the  structure  of  vortex  is  smeared  out.  In  the  new 
scheme,  on  the  other  hand,  the  strength  of  the  vortex  does  not  weaken  even  with  the  cell  length  being  ten  times  the  mean 
free  path.  At  Re- 200,  the  new  scheme  can  produce  the  vortex  shedding  behind  a  flat  plate  with  the  cell  dimension  of 
20  for  the  inclined  plate,  or  12^  in  the  case  of  a  normal  plate  to  the  flow  direction,  although  the  conventional  one 
cannot.  The  new  intermolecular  collision  scheme  is  verified  effective  for  unsteady  flows  such  as  vortex  shedding. 


INTRODUCTION 

The  DSMC  method  [1]  has  features  different  from  continuum  calculation  methods  like  the  finite  difference 
method.  The  DSMC  method  has  no  unstable  factor  that  a  calculation  diverges  and  has  the  merit  that  its  boundary 
condition  can  be  easily  determined.  Since  the  DSMC  is  a  microscopic  calculation  method  that  deals  with  molecules 
directly,  it  might  become  possible  to  elucidate  unstable  fluid  phenomena  such  as  a  turbulent  flow.  The  most  serious 
problem,  when  the  DSMC  method  is  applied  to  a  continuum  fluid,  is  that  it  is  necessary  to  divide  the  flowfield  into  a 
network  of  cells  the  dimension  of  which  is  smaller  than  the  local  mean  free  path.  At  the  meeting  of  the  RGD23,  one 
of  the  authors  has  proposed  a  new  intermolecular  collision  scheme  [2]  in  which  velocities  of  a  molecule  are 
modified  according  to  positions  of  a  collision  pair.  Although  the  scheme  was  verified  to  be  effective  with  a  one¬ 
dimensional  normal  shock  wave  and  with  a  two-dimensional  vortex  in  a  square  cavity,  the  original  version  of  the 
scheme  (Usami- system-0  or  U-sy stem-0)  was  not  perfect  for  an  axisymmetric  simulation  or  a  3-D  simulation  of  a 
supersonic  free  jet  and  the  velocity  modification  was  obliged  to  be  reduced  to  50%  of  the  full  value  because  the 
distortion  of  the  profile  appeared  with  the  full  modification.  Afterwards,  the  stability  of  the  calculation  was 
improved  (even  not  perfect)  with  the  addition  of  the  limitation  in  the  flow  velocity  ratio  and  the  temperature  ratio 
between  the  locations  of  two  molecules.  And  also,  it  was  found  that  the  use  of  the  bi-linear  or  the  tri-linear 
interpolation  for  calculation  of  the  temperature  and  the  flow  velocity  at  any  location  within  a  cell  can  reduce  the 
amount  of  the  storage  required  (U-system-1).  In  the  new  collision  scheme  the  total  momentum  or  the  total  energy  of 
each  collision  pair  is  not  necessarily  conserved  during  a  collision,  though  the  overall  value  of  all  molecules  can  be 
conserved  within  statistical  fluctuation.  After  various  researches  from  such  a  point  of  view,  we  have  found  the 
reliable  scheme  (U-sy stem-2)  in  which  the  stability  improves  much  by  an  additional  procedure,  and  the  full 
modification  of  the  velocity  becomes  possible  even  in  an  axisymmetric  simulation  or  a  3-D  simulation  of  a 
supersonic  free  jet  where  sudden  changes  occur  of  flow  properties.  The  new  collision  scheme  has  been  almost 
completed  for  steady  problems,  and  next,  it  is  necessary  for  us  to  investigate  an  effectiveness  of  the  scheme  to 
unsteady  problems.  We  target  a  generation  of  vortices  for  flow  past  a  flat  plate.  Under  low  Reynolds  number,  the 
steady  flow  past  a  plate  at  an  angle  of  attack  of  45  degrees  produces  one  vortex  behind  the  plate.  With  increasing  Re , 
however,  a  periodic  vortex  separation  occurs,  and  this  is  an  unsteady  state  problem.  We  also  apply  the  new  collision 
scheme  to  a  velocity  profile  in  a  circular  tube  under  the  condition  of  a  turbulent  flow  even  though  it  is  still  at  a 
testing  stage. 
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NEW  INTERMOLECULAR  COLLISION  SCHEME 


The  average  distance  between  molecules  of  a  collision  pair  is  in  proportion  to  the  cell  length  when  two  molecules 
are  chosen  at  random  among  molecules  in  the  cell.  If  the  cell  length  is  larger  than  the  local  mean  free  path,  the 
molecules  that  should  not  collide  normally  may  come  to  collide,  and  the  result  of  simulation  may  be  distorted.  It 
would  be  favorable  that  in  an  intermolecular  collision  the  second  molecule  as  a  collision  partner  could  be  selected  at 
the  same  location  that  the  first  molecule  is  selected  at.  In  the  new  scheme,  the  velocity  of  a  molecule  being  at  a 
different  location  is  modified  as  if  the  molecule  would  be  located  at  the  same  position  as  the  first  selected  molecule. 
We  assume  that  the  velocity  distribution  in  all  places  is  in  local  equilibrium,  and  we  call  two  molecules  P  and  Q  that 
participate  in  the  intermolecular  collision.  The  velocity  of  P  is  one  value  (point  A  in  Fig.  1)  of  the  velocity 
distribution  f\  at  the  location  of  P  and  the  velocity  of  Q  is  one  value  of  the  velocity  distribution  f2  at  the  location 
of  Q.  To  transfer  the  molecule  P  to  the  same  position  as  Q,  we  should  move  the  point  A  to  the  appropriate  position 
within  the  velocity  distribution  f2  .  The  most  natural  method  to  move  the  point  A  into  the  velocity  distribution  f2  is 
that  the  relative  location  of  P  within  f2  is  maintained  as  the  relative  location  A  within  fx  .  This  velocity 
modification  is  calculated  easily.  The  x  velocity  component  u2  in  f2  (the  point  B  in  Fig.l)  is  calculated  from  the 
original  velocity  ux  of  P  using  Eq.(l). 

u2=u2+{ul-ul)Jr^Jr~i  (i) 

where  Ux  and  TxX  are  the  flow  velocity  and  the  temperature  of  fx ,  and  U2  and  Tx2  are  those  of  f2  .  To  use  Eq.(l), 
it  is  necessary  to  find  the  temperature  and  the  flow  velocity  at  the  location  of  the  molecule  P  and  also  find  those  at 
the  location  of  the  molecule  Q.  To  obtain  the  temperature  and  the  flow  velocity  in  the  flowfield,  we  must  first 
calculate  the  flow  by  the  DSMC  method  using  the  conventional  collision  scheme  in  a  coarse  cell  network.  We  call 
this  calculation  Calcu_l.  The  variations  of  the  temperature  and  the  flow  velocity  (along  x,  y,  and  z  axes)  within 
the  cell  are  each  approximated  as  a  linear  function.  The  temperature  and  the  flow  velocity  at  any  location  are 
calculated  using  the  results  obtained  by  Calcu_l.  If  the  degree  of  nonequilibrium  is  small,  the  average  temperature 
may  be  used  such  as  Tx  =  Ty  =  TZ  =T  to  reduce  the  memory.  The  temperature  and  the  flow  velocity  obtained  by 
Calcu_l  must  be  corrected  by  the  calculation  that  follows  Calcu-1  and  so  on.  In  unsteady  problems,  especially,  the 
time  interval  between  Calcu_(i)  and  Calcu_(i+1)  should  be  small  enough  that  the  time  delay  would  not  have  a  bad 
influence  on  the  velocity  modification.  The  core  of  an  intermolecular  collision  calculation  is  the  same  as  the 
conventional  calculation,  although  using  the  modified  velocity  u2  (  v2,  w2  )  of  the  molecule  P  calculated  by  Eq.(l). 
As  a  result,  let's  assume  that  the  velocity  of  the  molecule  P  changes  to  the  point  C  (u2b)  on  the  velocity  distribution 
f2  as  in  Fig.l.  We  should  bring  the  post-collision  velocity  of  the  molecule  P  to  the  position  D  (  uXb  )  on  the  velocity 
distribution  fx  because  P  must  return  back  to  the  original  location.  The  velocity  modification  is  performed  by  Eq.(2) 
as  in  the  first  modification  where  Eq.(l)  is  used. 


Molecular  velocity  component 

FIGURE  1.  Modification  of  velocity  in  the  new 
intermolecular  collision  scheme 
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ulb  —  +  ( U2b  —U2  )y]Txi  /  Tx2  (2) 

Note  that  although  Eq.(l)  and  Eq.(2)  have  been  derived  with  the  assumption  of  equilibrium  distribution  functions, 
these  equations  themselves  and  the  new  collision  scheme  are  not  restricted  only  on  equilibrium  problems.  In  the 
third  version  of  the  new  collision  scheme  (U-system-2),  we  impose  the  limitation  of  the  ratio  between  flow 
velocities,  or  temperatures,  at  locations  of  two  molecules  as  mentioned  earlier.  We  also  use  a  bi-linear  or  a  tri-linear 
interpolation  to  estimate  flow  properties  at  any  location  within  a  cell.  For  example  of  the  bi-linear  interpolation  in  a 
two  dimensional  flowfield,  when  the  temperatures  of  four  corners  (x,  y),  (x+A,  y),  (x,  y+B )  and  (x+A,  y+B)  of  a 
rectangular  cell  are  7\,  r2,  T3  and  r4,  and  the  location  of  a  molecule  is  {x+a,  y+b ),  the  temperature  at  the  location  of 
the  molecule  is  calculated  by  Eq.(3). 

T=  (l-b/B){(l-a/A)  Tx  +  0 a/A )  T2 }  +  {b/B){{\- a/A)  T3+  ( a/A )  T4 }  (3) 

In  addition  to  that,  we  have  found  that  if  the  change  of  the  total  energy  during  a  collision  in  the  two  molecules 
exceeds  a  certain  limit,  for  example,  ±  10%  of  the  pre-collision  value,  the  probability  that  a  distortion  of  results  is 
caused  by  the  collision  pair  cannot  be  neglected  and  that  the  new  scheme  should  not  apply  to  the  intermolecular 
collision  of  such  a  pair.  Although  the  possibility  of  the  occurrence  of  the  collision  pair  is  very  weak,  for  example, 
less  than  0.0005,  and  an  addition  of  the  procedure  excluding  the  collision  pair  from  an  application  of  the  new 
scheme  does  not  affect  the  efficiency  of  the  new  scheme,  the  stability  improves  much  by  the  additional  procedure 
and  the  full  modification  of  the  velocity  becomes  possible  even  in  an  axisymmetric  simulation  or  a  3D  simulation  of 
a  supersonic  free  jet.  Figure  2  shows  the  comparison  of  density  profiles  of  the  axisymmetric  free  jet  (  Kn  —  3xlCT4 
and  pressure  ratio  50)  between  the  result  by  the  conventional  collision  scheme  with  249,600  cells  (or  15,600  cells) 
and  that  by  the  new  scheme  with  15,600  cells,  where  d  is  a  diameter  of  an  orifice.  The  new  scheme  realizes  an 
excellent  result  with  very  small  number  of  cells,  in  other  words,  with  very  small  number  of  molecules.  Anyone  who 
has  an  interest  in  the  new  collision  scheme  may  request  the  Fortran  source  code  (U-system-2)  of  the  free  jet 
calculation  to  the  author  through  E-mail  (usami@mach.mie-u.ac.jp).  Figure  3  shows  the  density  contours  in  two 
planes  normal  to  the  jet  axis  in  the  3-D  free  jet  analysis  ( Kn  —  3xlCT5  and  pressure  ratio  100)  calculated  with  the 
new  collision  scheme  and  the  conventional  one.  When  the  free  jet  flows  through  a  circular  orifice,  since  it  has 
instability  essentially,  the  axisymmetric  structure  is  not  maintained  and  a  petal  pattern  of  density  contours  is  known 
to  appear  in  a  plane  normal  to  the  jet  axis.  The  petal  patterns  in  the  jet  cross  sections  are  reproduced  very  clearly 
with  the  new  collision  scheme,  but  not  with  the  conventional  scheme,  where  13  xlO6  cells  and  36  xlO6  molecules 
are  used. 


-2-1012  z/d  -2-1012  z/d  -2-101  2  z/d  -2-101  2  z/d 

(a)  New  scheme  (b)  Conventional  scheme 


FIGURE  3.  Petal  patters  of  density  contours  in  planes  normal  to  the  jet  axis  (pressure  ratio  100) 


RESULT  AND  CONSIDERATION 


Steady  flow  past  an  inclined  flat  plate  with  one  vortex 

The  vortex  flow  behind  an  inclined  flat  plate  is  a  classical  problem  in  fluid  mechanics.  For  an  uncompressible 
fluid,  a  periodic  vortex  separation  begins  to  occur  in  the  vicinity  of  Reynolds  number  Re= 50.  Meiburg  adopted  the 
Re  number  in  his  paper  [3]  where  the  MD  method  and  the  DSMC  method  were  compared  in  respect  to  the 
manageability  of  vortices.  In  the  present  paper,  the  new  collision  scheme  and  the  conventional  scheme  of  the  DSMC 
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method  are  compared  with  respect  to  their  allowable  cell  dimension.  Figure  4  shows  the  two-dimensional  flowfield 
of  a  square  (580/1  X580>^ )  with  a  flat  plate  whose  length  is  60  A  (2.5  mm),  where  A  is  the  undisturbed  mean  free 
path.  A  flow,  the  speed  of  which  is  221  m/s  (Ma= 0.7),  enters  the  flowfield  at  an  angle  of  45  degrees  through  a  lower 
side  boundary  and  a  left-hand  side  boundary.  The  flowfield  is  divided  with  square  cells.  Three  types  of  linear 
dimension  of  a  cell  are  investigated:  A ,  5  A  and  10  A  .  The  VHS  molecule  of  argon  is  used.  The  pressure  is  0.9  Torr, 
the  temperature  is  288  K,  the  Kn  number  is  0.017,  and  the  Re  number  is  50.  The  total  number  of  molecules  is 
840,000  in  all  calculations.  Both  schemes  produce  a  similar  steady  vortex  behind  the  inclined  plate  when  the  cell 
length  is  A  (Fig.  5(a)).  In  the  conventional  collision  scheme,  however,  the  strength  of  vortex  becomes  weak  with  cell 
length  5  A  (Fig.  5(b)  upper  side)  and  the  vortex  disappears  with  cell  length  10  A  (Fig.  5(c)  upper  side).  On  the  other 
hand,  in  the  new  scheme,  the  same  strength  of  vortex  as  in  the  case  of  cell  length  A  is  maintained  even  when  the 
cell  length  becomes  10  A  (Fig.  5  lower  side). 


FIGURE  4.  Two  dimensional 
flowfield  used  to  simulate 
vortex  flows 


(a)  Cell  length  =  X  (b)  Cell  length  =  5  X  (c)  Cell  length  =  10  X 


FIGURE  5.  Steady  vortex  behind  an  inclined  flat  plate 


FIGURE  6.  Unsteady  vortex  flow  past  an  inclined  flat  plate 


FIGURE  7.  Vorticity  contours  behind 
an  inclined  plate  by  the  new  scheme 


Unsteady  flow  past  an  inclined  flat  plate  with  vortex  shedding 

The  steady  vortex  behind  the  inclined  flat  plate  has  been  obtained  at  Re= 50  as  mentioned  above.  To  produce  an 
unsteady  vortex  shedding.  Re  number  should  be  increased.  For  Re= 200  (Kn= 0.0043),  the  two  schemes  are  compared 
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in  the  same  flowfield  as  above-mentioned.  The  flowfield  is  divided  with  13,456  (116X  116)  square  cells.  The  linear 
dimension  of  a  cell  is  20  X  and  the  length  of  the  plate  is  240  A, .  The  pressure  is  3.5  Torr.  The  number  of  molecules 
is  840,000.  In  the  new  scheme,  the  flow  properties  (flow  velocity  and  temperature)  are  needed  to  modify  molecular 
velocities.  Since  they  must  change  with  time  in  case  of  unsteady  flows,  samplings  of  flow  properties  should  be  done 
at  a  sufficiently  short  time  interval  compared  with  the  time  scale  for  changing  speed  of  the  flowfield.  The  present 
sampling  interval  is  selected  to  be  about  one-twentieth  of  the  characteristic  time  for  a  periodic  vortex  separation. 
Figure  6  shows  the  results  of  streamline  obtained  with  the  different  collision  schemes.  The  new  collision  scheme  can 
produce  a  distinct  vortex  shedding,  although  the  conventional  scheme  cannot.  Figure  7  indicates  vorticity  contours 
obtained  by  the  new  collision  scheme.  A  negative  vorticity  (thin  lines)  is  produced  at  the  left-hand  side  of  the  plate 
and  a  positive  vorticity  (thick  lines)  is  produced  at  the  right-hand  side  periodically,  and  they  flow  to  the  downstream 
with  the  passage  of  time.  The  calculation  time  of  the  new  scheme  in  this  case  is  about  40%  greater  than  that 
compared  to  the  conventional  scheme.  It  has  the  inclination  to  increase  with  increasing  density. 

Vortex  shedding  behind  a  flat  plate  perpendicular  to  the  flow 

Periodic  vortex  shedding  behind  a  flat  plate  perpendicular  to  the  flow  is  also  investigated  with  the  two  collision 
schemes.  Figure  8  shows  the  rectangular  flowfield  that  is  divided  with  20,000  square  cells  (200 X  100).  The  linear 
dimension  of  the  cell  is  12  A  ,  and  the  length  of  the  flat  plate  is  240  A  .  The  number  of  molecules  is  640,000.  The 
pressure  of  argon  gas  is  4.3  Torr,  the  temperature  is  288  K,  the  flow  velocity  is  221  m/s  (Ma= 0.7),  Re= 200,  and 
Kn= 0.0043.  A  monitoring  cell  in  which  an  oscillation  of  flow  properties  is  investigated  is  located  at  the  center  of  the 
flowfield.  Although  only  small  flow  disturbances  behind  the  plate  are  found  in  the  result  of  the  conventional 
calculation  (Fig.  9),  regular  periodic  vortex  shedding  is  clearly  seen  as  in  Fig.  10  when  using  the  new  scheme.  The 
clockwise  vortex  (negative  vorticity)  is  produced  at  an  upper  edge  of  the  plate  and  the  counter  clockwise  vortex 
(positive  vorticity)  is  produced  at  a  lower  edge,  alternately.  They  flow  to  the  downstream  keeping  a  certain  distance 
between  them.  Figure  11  shows  vorticity  contours  by  the  new  scheme  (thin  lines  are  negative  and  thick  lines  are 
positive  vorticity).  An  oscillation  frequency  of  the  vortex  shedding  can  be  obtained  through  investigating  temporal 
variation  of  flow  properties,  especially  using  velocity  component  normal  to  the  flow  direction,  in  the  monitoring  cell. 
Strouhal  number  St  turns  out  to  be  about  0.2,  which  is  a  reasonable  value  according  to  the  published  data. 


FIGURE  8.  Flowfield  used  to  simulate  vortex 
shedding  behind  a  normal  flat  plate 


FIGURE  9.  Small  flow  disturbances  behind  a  normal  flat  plate  by 
the  conventional  collision  scheme 


FIGURE  10.  Regular  vortex  shedding  behind  a  normal  flat 
plate  by  the  new  collision  scheme 


FIGURE  11.  Vorticity  contours  behind  a 
normal  flat  plate  by  the  new  collision  scheme 
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Velocity  profile  in  a  circular  tube  under  a  turbulent  flow  condition 


Molecular  simulation  on  transition  of  flow  from  laminar  to  turbulent  is 
one  of  the  most  interesting  problems.  A  deformation  from  the  velocity 
profile  of  Poiseuille  flow  in  a  cylindrical  tube  is  investigated  with  the  new 
collision  scheme.  This  attempt  is  not  necessarily  performed  by  perfect 
algorithm  based  on  the  theoretical  endorsement  and  it  is  still  at  a  testing 
stage.  A  cross-section  (y-z  plane)  of  a  circular  tube  is  enclosed  by  a 
circumscribed  quadrate  that  is  divided  by  160,000  (400X400)  square  cells. 

Only  one  cell  is  available  in  the  v  direction  (the  flow  direction)  and  the 
periodic  boundary  condition  is  applied  to  the  boundary  perpendicular  to  the 
v-axis.  The  number  of  molecules  of  argon  gas  is  about  1260,000.  A  diameter 
of  the  tube  is  4  mm  and  an  average  flow  speed  is  69  m/s  (0.2  times  the  most 
probable  molecular  speed).  An  initial  velocity  profile  is  parabolic.  A  flow 
speed  would  decrease  with  time  because  the  flow  suffers  pipe  friction.  In  this 
calculation,  then,  a  momentum  loss  parallel  to  the  pipe  wall  is  compensated 
by  means  of  giving  small  momentum  to  all  molecules  equally  to  conserve 
the  total  value.  Moreover,  the  translational  energy  would  be  remarkably 
accumulated  with  increasing  time,  especially  in  the  new  collision  method, 
and  the  temperature  would  change  along  the  radial  direction  in  the  tube 
considerably.  To  prevent  this  undesirable  effect,  the  total  energy  over  all 
molecules  is  counted  just  after  the  intermolecular  collision  routine  and  a 
small  energy  is  subtracted  from,  or  against  long  odds  added  to,  the  energy  of 
each  molecule  equally  every  time  step,  so  that  the  total  energy  can  be 

conserved.  Enough  consideration  has  not  yet  been  done  about  an  adverse  effect  that  these  artificial  operations  cause 
for  this  simulation.  The  time  interval  between  updates  of  flow  properties  that  are  used  for  modification  of  molecular 
velocities  is  3000  times  the  mean  free  time  of  molecular  collisions.  Figure  12  shows  velocity  profiles  in  the  case  of 
pressure  400  Torr  ( Re=  10,000).  Though  the  departure  of  the  velocity  from  the  curve  of  laminar  flow  is  apparent, 
since  the  present  result  is  sensitive  to  the  time  interval,  more  careful  searches  about  it  will  be  needed  in  the  future.  In 
the  conventional  scheme,  on  the  other  hand,  all  profiles  of  flow  velocity  become  parabolic  irrespective  of  Re. 


FIGURE  12.  Velocity  profile 
in  a  circular  tube 


CONCLUSION 

The  present  paper  is  related  to  the  new  intermolecular  collision  scheme  of  the  DSMC  method  and  the  following 
conclusions  have  been  obtained. 

(1)  It  is  useful  for  an  improvement  of  stability  in  the  new  collision  scheme  that  the  new  scheme  applies  only  to 
the  collision  pair  whose  post-collision  energy  is  within  90%  to  110%  of  the  pre-collision  value. 

(2)  Although  the  strength  of  steady  vortex  behind  an  inclined  plate  becomes  weak  with  cell  length  5  A  and  the 
vortex  disappears  with  cell  length  10  A  in  the  conventional  scheme  at  Re= 50,  the  same  strength  of  vortex 
as  in  the  case  of  cell  length  A  is  maintained  in  the  new  scheme  even  when  the  cell  length  becomes  10  A  . 

(3)  In  the  case  of  cell  length  20  A  ,  the  new  collision  scheme  can  produce  a  distinct  vortex  shedding  behind  an 
inclined  flat  plate  at  Re— 200,  although  the  conventional  scheme  cannot. 

(4)  Regular  periodic  vortex  shedding  is  clearly  seen  behind  a  flat  plate  normal  to  the  flow  direction  and 
Strouhal  number  about  0.2  is  obtained  at  Re= 200  when  using  the  new  scheme  with  cell  length  12  A  , 
although  only  small  flow  disturbances  are  found  in  the  result  of  the  conventional  scheme. 
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